8:7 Journal of Turbulence eigenPipeArc 



Journal of Turbulence 

Vol. , No. , Jan 2007, 1-28 



Dynamical Eigenfunction Decomposition of 
Turbulent Pipe Flow 

Andrew Duggleby,^ Kenneth S. Ball,^ Mark R. Paul,t and Paul F. Fischer§ 

f Department of Mechanical Engineering, Virginia Polytechnic Institute and State 
University, Blacksburg, VA 24061 
§ Mathematics and Computer Science Division, Argonne National Laboratory, Argonne, IL 

60439 



(in review) 



The results of an analysis of turbulent pipe flow based on a Karhunen-Loeve decomposition are 
presented. The turbulent flow is generated by a direct numerical simulation of the Navier-Stokes 
equations using a spectral element algorithm at a Reynolds number Re T = 150. This simulation 
yields a set of basis functions that captures 90% of the energy after 2,453 modes. The eigenfunctions 
are categorised into two classes and six subclasses based on their wavenumber and coherent vorticity 
structure. Of the total energy, 81% is in the propagating class, characterised by constant phase 
speeds; the remaining energy is found in the non propagating subclasses, the shear and roll modes. 
The four subclasses of the propagating modes are the wall, lift, asymmetric, and ring modes. The wall 
modes display coherent vorticity structures near the wall, the lift modes display coherent vorticity 
structures that lift away from the wall, the asymmetric modes break the symmetry about the axis, 
and the ring modes display rings of coherent vorticity. Together, the propagating modes form a wave 
packet, as found from a circular normal speed locus. The energy transfer mechanism in the flow is 
a four step process. The process begins with energy being transferred from mean flow to the shear 
modes, then to the roll modes. Energy is then transfer ed from the roll modes to the wall modes, 
and then eventually to the lift modes. The ring and asymmetric modes act as catalysts that aid in 
this four step energy transfer. Physically, this mechanism shows how the energy in the flow starts at 
the wall and then propagates into the outer layer. 

Keywords: Direct numerical simulation, Karhunen-Loeve decomposition, turbulence, pipe flow, 
mechanism 



1 Introduction 

Turbulence, hailed as one of the last major unsolved problems of classical 
physics, has been the subject of numerous publications as researchers seek to 
understand the underlying physics, structures, and mechanisms inherent to the 
flow [1]. The standard test problem for wall-bounded studies historically has 
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Table 1. Summary of existing algorithms for the DNS of turbulent pipe flow, 
citing the first reported use for each algorithm. As seen, most approaches use a 
2nd-order radial discretization. The reason is that the standard spectral meth- 
ods, in the presence of the coordinate singularity, achieve only 2nd-order con- 
vergence instead of geometric convergence. 

Authors Radial Discretization 

Eggels et al. (1994) [8] 2nd-order FD - first pipe DNS 

Verzicco and Orlandi (1996) [9] 2nd-order FD flux based 

Fukagata and Kasagi (2000) [10] 2nd order FD - energy conservative 

Shan et al. (1999) [6] "Piecemeal" Chebyshev 

Loulou et al.(1997) [7] B-splinc 



been turbulent channel flow because of its simple geometry and computational 
efficiency. Even though much insight has been achieved through the study of 
turbulent channel flow, it remains an academic problem because of its infinite 
(computationally periodic) spanwise direction. 

The next simplest geometry is turbulent pipe flow, which is of interest be- 
cause of its real-world applications and its slightly different dynamics. The 
three major differences between turbulent pipe and channel flow are that pipe 
flow displays a log layer, but overshoots the theoretical profile until a much 
higher Reynolds numbers (Re = 3000 for a channel versus Re = 7442 for 
a pipe [2, 3]), has an observed higher critical Reynolds number, and is lin- 
early stable to an infinitesimal disturbance [4, 5] . Unfortunately, few direct 
numerical simulations of turbulent pipe flow have been carried out because 
of the complexity in handling the numerical radial singularity at the origin. 
Although the singularity itself is avoidable, its presence causes standard high- 
order spectral methods to fail to converge exponentially. As a result, only a 
handful of algorithms found in the literature for turbulent pipe flow; the first 
reported use for each algorithm is listed in Table [TJ These algorithms typically 
use a low-order expansion in the radial direction. Only Shan et al. [6] using 
concentric Chebyshev domains ( "piecemeal" ) and Loulou et al. [7] using basis 
spline (B-spline) polynomials provide a higher-order examination of turbulent 
pipe flow. Using a spectral element method, this study provides not only a 
high-order examination but the first exponentially convergent investigation of 
turbulent pipe flow through direct numerical simulation (DNS). 

With this DNS result, one of the studies that can be performed with the full 
flow field and time history it provides is an analysis based on an orthogonal 
decomposition method. In such methods, the flow is expanded in terms of 
a natural or preferred turbulent basis. One method used frequently in the 
field of turbulence is Karhunen-Loeve (KL) decomposition, which uses a two- 
point spatial correlation tensor to generate the eigenfunctions of the flow. 
This is sometimes referred to as proper orthogonal decomposition, empirical 
orthogonal function, or empirical eigenfunction analysis. 

Work in this area was pioneered by Lumley, who was the first to use the KL 
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method in homogeneous turbulence [11, 12]. This method was later applied to 
turbulent channel flow in a series of papers by Ball, Sirovich, and Keefe [13, 14] 
and Sirovich, Ball, and Handler [15], who discovered plane waves and propa- 
gating structures that play an essential role in the production of turbulence 
through bursting or sweeping events. To study the interactions of the propagat- 
ing structures, researchers have examined minimal expansions of a turbulent 
flow [16-19]. These efforts have led to recent work by Webber et al. [20,21], 
who examined the energy dynamics between KL modes and discovered the 
energy transfer path from the applied pressure gradient to the flow through 
triad interactions of KL modes. 

This present study uses a spectral element Navier-Stokes solver to generate a 
globally high-order turbulent pipe flow data set. The Karhunen-Loeve method 
is used to examine the turbulent flow structures and dynamics of turbulent 
pipe flow. 



2 Methodology 

The direct numerical simulation of the three-dimensional time dependent 
Navier-Stokes equations is a computationally intensive task. By fully resolving 
the necessary time and spatial scales of turbulent flow, however, no subgrid 
dissipation model is needed, and thus a turbulent flow is calculated directly 
from the equations. DNS has one main advantage over experiments, in that 
the whole flow field and time history are known, enabling analyses such as the 
Karhunen-Loeve decomposition. 

Because of the long time integration and the grid resolution necessary for 
DNS, a high-order (typically spectral) method is often used to keep numerical 
round-off and dissipation error small. Spectral methods and spectral elements 
use trial functions that are infinitely and analytically differentiable to span the 
element. This approach decreases the global error exponentially with respect 
to resolution, in contrast to an algebraic decrease with standard methods such 
as finite difference or finite element methods [22]. 

2.1 Direct Numerical Simulation 

2.1.1 Numerical Methods. This study uses a spectral element Navier- 
Stokes solver that has been developed over the past 20 years [23, 24] to solve 
the Navier-Stokes equations: 



d t U + U • VU = -VP + Re-^U 



(1) 
(2) 



V • U = 0. 
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Figure 1. Cross section of the spectral element grid for pipe flow. The spectral element algorithm 
avoids the singularity found in polar cylindrical coordinates. 

In the above equations, U = (U r , Uq, U z ) is the velocity vector corresponding to 
the radial, azimuthal, and streamwise direction respectively; Re T = U T R/v is 
the Reynolds number; R is the radius of the pipe; v is the kinematic viscosity; 
and U T = \/r w /p is the shear velocity based upon the wall shear stress t w 
and density p. This solver employs a geometrically flexible yet exponentially 
convergent spectral element discretization in space, in this case allowing the 
geometry to be fitted to a cylinder. The domain is subdivided into elements, 
each containing high-order (usually 11-13) Legendre Lagrangian interpolants. 
The resulting data localisation allows for minimal communication between 
elements, and therefore efficient parallelization. Time discretization is done 
with third-order operator splitting methods, and the remaining tensor-product 
polynomial bases are solved by using conjugate gradient iteration with scalable 
Jacobi and hybrid Schwarz/multigrid preconditioning [25]. 

2.1.2 Geometry. Spectral elements are effective in cylindrical geometries 
[26] and their use elegantly avoids the radial singularity at the origin. The 
mesh is structured as a box near the origin and transitions to a circle near 
the pipe walls (Figured]), maintaining a globally high-order method at both 
the wall and the origin. In addition to avoiding the numerical error associated 
with the singularity, the method also avoids the time-step restriction due to 
the smaller element width at the origin of a polar-cylindrical coordinate sys- 
tem, which could lead to potential violations of the Courant-Friedrichs-Levy 
stability criteria. 

Each slice, as shown in Figure [IJ has 64 elements, and there are 40 slices 
stacked in the streamwise direction, adding up to a length of 20R. Each element 
has 12th-order Legendre polynomials in each direction for Re r = 150. This 
discretization results in 4.4 million degrees of freedom. Near the wall, the 
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grid spacing normalised by wall units (v/U T ) denoted by the superscript + is 
Ar + ~ 0.78 and (RA9) + ~ 4.9, where r is the radius and 6 is the azimuthal 
angle. Near the centre of the pipe, the spacing in Cartesian coordinates is 
A + ~ 3.1. The streamwise grid spacing is a constant Az + = 6.25 throughout 
the domain where z is the streamwise coordinate. 



2.1.3 Benchmarking at Re x = 180. Benchmarking was performed at 
Re T = 180 with the experiments and DNS of Eggels et al. [8] and the DNS of 
Fukagata and Kasagi [10]. For this higher Reynolds number flow, 14th-order 
polynomials were used, giving grid spacings near the wall of Ar + ~ 0.80 and 
(RA6) + sa 5.0, A+ w 3.2 in the centre, and Az + = 6.42. Eggels et al. and 
Fukagata and Kasagi used a spectral Fourier discretization in the azimuthal 
and axial directions and then a 2nd-order finite difference discretization in the 
radial direction. Also, both groups used a domain length of 10R and grid sizes 
in their DNS studies of 96 x 128 x 256 for r, 6, and z directions, respectively. 

The turbulent flow was tripped with a solenoidal disturbance, and the simu- 
lation was run until transition was complete. This was estimated by examining 
the mean flow and root-mean-squared statistics over sequential periods of 100 
t + until a statistically steady state was achieved, marked by a change of less 
than 1 percent between the current and previous period. In this case, it took 
7000i + to arrive at a statistically steady state. Convergence was checked by 
examining the total fluctuating energy after 1000 time steps for increasing 
polynomial orders 8, 10, 12, and 14 with a fixed number of elements (2560). 
The absolute value of the difference of the total fluctuating energy for each 
simulation and the value obtained using 16th order polynomials was plotted 
versus the total number of degrees of freedom in Figure [2J The exponential 
decay of the error with increasing degrees of freedom shows that our spectral 
element algorithm achieves geometric or exponential convergence. 

The mean flow profiles in Figure [3] correspond well with the hot wire 
anemometry (HWA) results of Eggels et al. [8], but as seen in the root-mean- 
squared (rms) statistics shown in Figure [H the spectral element calculation 
shows a lower peak U z ^ rms and higher peak Ue, r ms and U r)rms compared to 
the HWA, particle image velocimetry (PIV), and laser Doppler anemometry 
(LDA) results. These results are in contrast to those with channel flow, as 
reported by Gullbrand [27], where the 2nd order finite difference methods 
undershoot the spectral method wall-normal velocity rms. 

When compared to the experimental results of Eggels et al. [8], the spectral 
element results are in better agreement than the 2nd-order finite difference 
results. We note that Eggels et al. report that their (PIV) system had diffi- 
culties capturing U r near the wall and near the centerline of the pipe due to 
reflection, which could explain the deviation of all of the DNS results in that 
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10° 




Degrees of freedom 

Figure 2. Convergence plot showing the error in fluctuating total energy of the flow versus degrees 
of freedom after 1000 time steps for a fixed grid and polynomial orders of 8, 10, 12, and 14 (circles), 
using 16th order polynomials to approximate the exact solution. The exponentially decay of the 
curve fit (solid) shows that our spectral element algorithm is geometrically or exponential 

convergent. 

area. 

A second major difference our spectral element method and the previous 
pipe DNS using 2nd-order finite difference is the domain size. With the spec- 
tral element method, which results in a global high-order convergence, a do- 
main size of 10R yielded unphysical results in the flow, as seen in the bulge in 
the azimuthal Ug rms in Figure [5] at y + = 55 and y + = 120. This unphysical 
result arose even with a more refined grid, and is therefore not a function of 
under-resolution. However, when the domain of the spectral element method 
was extended to 20R, this problem disappeared. We surmise that the 2nd- 
order finite difference method dissipated some large-scale structures after 10R 
that the higher-order spectral element case appropriately resolves. This undis- 
sipated structure, because of the periodic boundary conditions, then re-enters 
the inlet and causes the unphysical bulges in the azimuthal rms profile. This 
result is also supported by the work of Jimenez [28] in turbulent channel flow 
(Re r =180) that shows large-scale structures that extend well past the domain 
size of z + = 1800 = WR. 

This benchmark confirms that the spectral element algorithm, at the given 
grid resolution and domain size, will generate the appropriate turbulent flow 
field and time history to perform the KL decomposition. 

2.2 Karhunen-Loeve Decomposition 



For completeness, the Karhunen-Loeve Decomposition method is briefly de- 
scribed here, but for more detail, see [11-16,29-32]. 



February 2, 2008 



i:7 Journal of Turbulence eigenPipeArc 



A. Duggleby, K.S. Ball, M.R. Paul, and P.F. Fischer 




Figure 3. Benchmark of mean velocity profile for the spectral element algorithm vs low-order 
methods of Eggels et al. [8] and Fukagata and Kasagi [10] for Re T = 180. The theoretical line is the 
law of the wall U + = y + and the log layer U + = 1/0.41 log y + + 5.5. Deviations from the log layer 
are expected in turbulent pipe flow until much higher Reynolds numbers. 

2.2.1 Numerical Method. The Karhunen-Loeve method is the solution of 
the two-point velocity correlation kernel equation defined by 



/ / / K(x,x')$>(x.')r'dr'd9'dz' = A3> 
Jo Jo Jo 

K(x,x') = (u(x)®u(x / )), 



(3) 
(4) 



where u(x) = U(x) — U(x) is the fluctuating component of the velocity and 
where the mean U(x) is determined by averaging over both homogeneous 
planes and time. The angle brackets represent an average over many time steps, 
on the order of t + ~ 40, 000 ~ 50, 000, to sample the the entire attractor. The 
(g> denotes the outer product, establishing the kernel as the velocity two-point 
correlation between every spatial point x = (r, 9, z) and x' = (r', 9', z'). 

For turbulent pipe flow, with two homogeneous directions providing trans- 
lational invariance in the 9 (azimuthal) and z (streamwise) direction, Eq. ^) 
becomes 



K(x, x) = K(r, r', 9-9',z- z') 



(5) 
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Figure 4. Benchmark of the rms profiles for Re T =180 comparing the spectral element algorithm vs 
the low-order methods of Eggels et al. [8] and Fukagata and Kasagi [10]. As seen near y+ = 20 in 
the w^ rms and near jj+ = 50 in the u^ rms , the higher-order method is closer to the experimental 

results than the 2nd-order methods. 

= lC{n,m;r,r')e in8 e i27rmz/L . (6) 

Thus, given the kernel in Eq. ([5]), the eigenfunctions have the form 

* (m>n) (r,9,z) = *(n,m;r)e m V 2 ™ z / L , (7) 

where n is the azimuthal wavenumber and m the streamwise wavenumber. 
The determination of \I/ is then given by 

f R 

/ K(m,n;r, r')*&*(m, n;r r )r'dr = X mn ^(m,n;r), (8) 
Jo 

/C(m, n; r, r ) = (u(n, m; r) ® u(n, m; r )) (9) 

where the * denotes the complex conjugate, A is the eigenvalue, and u(n, m; r) 
is the Fourier transform of the fluctuating velocity in the azimuthal and axial 
direction. 

In the present problem, 2,100 snapshots of the flow field were taken, cor- 
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Figure 5. Effect of a short domain (10R) with a high-order method, showing unphysical deviations 
in the rms profiles near y+ = 50 and y + = 120 for the spectral element algorithm for Re T = 180. A 
bulge, highlighted by the red ovals, is present in the azimuthal rms that is not found in the 
low-order methods of Eggels et al. [8] and Fukagata and Kasagi [10] . By extending the domain to 

20R, the unphysical bulge is removed, bflt is therefore surmised that the unphysical result is 
possibly due to the more dissipative nature of the low-order schemes that washes away a long-scale 
structure that the present high-order scheme captures and re-enters the inlet from the periodic 

outflow conditions. 



responding to one snapshot every eight viscous time steps ( t + = U^t/u). 
The results of each snapshot were projected to an evenly spaced grid with 
101 x 64 x 400 points in r, 9, and z, respectively. The Fourier transform of 
the data was then taken and the kernel assembled. This kernel was averaged 
over every snapshot to generate the final kernel to be decomposed. Since the 
dimension of K is 303 (given by three velocity components on 101 radial grid 
points) there are 303 eigenfunctions and eigenvalues for each Fourier wavenum- 
ber pair (m, n). The eigenfunctions are ranked in descending order with the 
quantum number q to specify the particular eigenfunction associated with the 
corresponding eigenvalue. Thus, it requires a triplet (m, n, q) to specify a given 
eigenfunction. 

The eigenfunctions ty q (m,n;r) are complex vector functions and are nor- 
malised so that the inner product of the full eigenfunction 3?( mjn (J ) is of unit 
length, namely (&( m ,n,q)>®(n',m',q')) = <W/<W<W, where 5 is the Kronecker 
delta. The eigenvalues physically represent the average energy of the flow cap- 
tured by the eigenfunction 3?( m . nj(? ), 
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\m,n,q) = (\(™, $ (m,n,g))\ 2 ) ■ (10) 

We note that the eigenfunctions, as an orthogonal expansion of the flow field, 
retain the properties of the flow field, such as incompressibility and boundary 
conditions of no slip at the wall. 

2.2.2 Symmetry Considerations. The pipe flow is statistically invariant 
under azimuthal reflection, 

Re : (r, 9, z, u r , u e ,u z ) -> (r, -9, z, u r , -u e ,u z ), (11) 

and taking advantage of it reduces the total number of calculations as well as 
the memory and storage requirements. We note that, because of its geometry, 
turbulent channel flow has two more symmetries - a vertical reflection, and a 
x-axis rotation - that are not present in the pipe, since a negative radius is 
equivalent to a 180 degree rotation in the azimuthal direction. 

A major consequence of this symmetry is that the resulting eigenfunctions 
are also symmetric, and the modes with azimuthal wave number n will be the 
azimuthal reflection of the modes with wave number — n, 

■ {^(m,n,q)'^(m,n,g)i^(m,n,q)) ~^ (^(m,—n,q)i ~ ^(m,—n,q)> ^(m,—n,q))i (1^) 

thus reducing the total computational memory needed for this calculation. 

2.2.3 Time- Dependent Eigenfunction Flow Field Expansion. The KL 

method provides an orthogonal set of eigenfunctions that span the flow field. 
As such, the method allows the flow field to be represented as an expansion 
in that basis, 

u(r,9,z,t)= Q(m,n,q) (0* ), (13) 

(m,n,q) 

with 

"'(m,n,g) 
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Since the Fourier modes are orthogonal to each other, equation 1141 becomes 



with u being the Fourier transform of u in the azimuthal and streamwise 
direction with wavenumbers n and m, respectively. 

The time history of the eigenfunctions can be used to examine their inter- 
actions, such as the energy interaction examined by Webber et al. [21] and 
bursting events by Sirovich et al. [15]. 

3 Results 

This section presents the results of our analysis of turbulent pipe flow based 
on KL decomposition. 

3.1 Mean Properties of Flow and Flow Statistics 

DNS was performed for Re T = 150. The Reynolds number based on the cen- 
treline velocity is Re c ~ 5700 and based on the mean velocity is Re m ~ 4300. 
This range is above the observed critical Reynolds number for pipe flow and 
exhibits self-sustaining turbulence, as seen from the fluctuating time history 
of the mean velocity shown in Figure El The velocity profile, seen in Figure 
shows the mean velocity with respect to wall units away from the wall 
{y + = {r — R)U T /u). The profile fits the law of the wall but fails to conform 
to the log law, and as mentioned in Section 1, the log law is not expected for 
turbulent pipe flow until much higher Reynolds number. 

The rms velocity fluctuation profiles and the Reynolds stress profile is shown 
in Figures [Ha] and I8b[ The streamwise fluctuations peak near y + = 16. The 
azimuthal and axial velocities show a weaker peak near y + = 39 and y + = 55, 
respectively, and then remain fairly flat throughout the pipe. The Reynolds 
stress v z v r has a maximum of 0.68 at y + = 31. 

3.2 Eigenvalue Spectrum 

As discussed in Section 2.2, the eigenvalues represent the energy of each eigen- 
function. By ordering the eigenvalues from largest to smallest, one can min- 
imise the number of eigenfunctions, N, needed to capture a given percentage 
of the energy of the flow. Table [2] shows the first 25 eigenfunctions, and Figure 
[9]shows the running total of energy versus modes. The first mode to have q ^ 1 




(15) 
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Figure 6. Variations in bulk mean velocity u m (t + ) consistent with turbulent pipe flow. The 
dashed line is the time averaged mean velocity. 




Figure 7. Mean velocity profile (solid) compared to the theoretical law of the wall u + = y + and 
the log layer u + = 1/0.41 log y + + 5.5 (dashed). As in the Re T = 180 case, the overshoot of the log 
layer above the theoretical profile is expected in turbulent pipe flow until much higher Reynolds 

number. 
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Figure 8a. Root-mean-square statistics across the 
pipe for the radial (solid), azimuthal (dashed), and 
streamwise (dot-dashed) velocities. The results show 
the standard profiles, similar to the Re T = 180 rms 
profiles. 



Figure 8b. Reynolds stress v z v r across the pipe 
(solid), peaking at 0.68 at y+ = 31, and the total 
shear stress distribution (dashed). 




2000 3000 
Eigenfunction 



5000 



Figure 9. Percentage of energy retained in the KL expansion with the 90% mark achieved after 
2453 modes (Karhuncn-Loeve dimension). 



is the 41st mode (0,0,2). The 90% mark is reached at D KL = 2,453, where 
Dkl can be considered as a measure of the intrinsic dimension of the chaotic 
attractor describing the turbulence as discussed by Zhou and Sirovich [17,32]. 
These eigenfunctions are the preferred or natural basis function for turbulent 
pipe flow, and insight is gained by observing their qualitative structure. 
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Tabic 2. First 25 eigenvalues ranked in descending order of 
energy. 



Index 






c l 


H 1 nV^Tl "\7"£* 1 1 1 O 
J_J1H,L11 VCLl Lit; 


P-npro-v f% Totals 


1 





6 


1 


1.61 


2.42 % 


2 





5 


1 


1.48 


2.22 % 


3 





3 


1 


1.45 


2.17 % 


4 





1 


1 


1.29 


1.93 % 


5 





2 


1 


1.26 


1.88 % 


6 


1 


5 


1 


0.936 


1.40 % 


7 


1 


6 


1 


0.917 


1.37 % 


8 


1 


3 


1 


0.902 


1.35 % 


9 


1 


1 


1 


0.822 


1.23 % 


10 





1 


1 


0.805 


1.20 % 


11 


1 


7 


1 


0.763 


1.14 % 


12 


1 


2 


1 


0.683 


1.02 % 


13 





7 


1 


0.646 


0.97 % 


14 


2 


1 


1 


0.618 


0.92 % 


15 





8 


1 


0.601 


0.90 % 


16 


2 


5 


1 


0.580 


0.87 % 


17 


1 


1 


1 


0.567 


0.85 % 


18 


2 


7 


1 


0.524 


0.78 % 


19 


1 


8 


1 


0.483 


0.72 % 


20 


2 


6 


1 


0.476 


0.71 % 


21 


2 


.3 


1 


0.454 


0.68 % 


22 


2 


2 


1 


0.421 


0.63 % 


2.3 


2 


8 


1 


0.375 


0.56 % 


21 


1 


9 


1 


0.358 


0.54 % 


25 


3 


1 


1 


0.354 


0.53 % 



3.3 Structure of the Eigenfunctions 

The eigenfunctions resulting from the KL decomposition can be categorised 
into two distinct classes and six sub-classes based on their characteristics, 
as listed in Table [3] and depicted in Figure QJH It is observed that certain 
characteristics of the eigenfunction are dependent on the wavenumber. The 
modes with higher azimuthal than axial wavenumber turn more than they lift, 
and as such the near-wall stretching structures are found. Likewise the modes 
with higher axial than azimuthal wavenumber lift more than they turn, and as 
such the lifting structures that extend from the near wall to the outer region 
are found. 

This eigenfunction expansion, using both their structure and their time- 
dependent coefficients, provides a framework for further analysis and com- 
parison and has led to understanding the mechanism of drag reduction by 
spanwise wall oscillation [36] and the mechanism of relaminarization [37]. 

In order to show the qualitative and quantitative classifications of these 
eigenfunctions, the most energetic function of each class is presented in Figures 
lllal through I16dl For each class, subfigures (a) show iso-contours and slices 
of coherent vorticity, and a slice displaying the Reynolds stress and coherent 
vorticity together. The coherent vorticity is defined as the rotating part of the 
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Tabic 3. Structures of turbulent pipe flow as classified by wavenumbcr. 




Struct urc 


13efi tilt io n 


Energy 


Descript ion 


Figure 


Propagating Modes 


(m,n,q), 


80.58% 


Constant phase speed 




(a) Wall 


n > m 


35.22% 
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Near wall streamwise vortices 


Fig.lTSalf 
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1.08% 


Non-zero centerline streamwise velocity 


Fig.[T6a}d 








No coherent vorticity 






Figure 10. Subclass classification according to streamwise (m) and spanwise (n) wavenumbcrs. 



fluid, defined as the imaginary component of the eigenvalues of the velocity 
strain rate tensor, following the work of Chong et al. [33]. The next four 
subfigures for each class are the plots of the eigenfunctions, both the real 
(b) and imaginary component (d), and the time history amplitude (c) and 
phase (e). To complete the class description, the coherent vorticity of two 
more representative eigenfunctions are shown in subfigures (f) and (g). For 
the shear modes (Figure [16a]) , since they do not have coherent vorticity, only 
the eigenfunction and time history are plotted. 
The first subclass of interest is the wall eigenmodes. The coherent vorticity 
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plots of three example wall modes are seen in Figures lllal lllfl and |llg| for 



the (1,5,1), (1,3,1), and (2,4,1) modes, respectively. Each mode consists of a 
travelling-wave coherent vortex near the wall, and all eigenmodes with greater 
spanwise wavenumber than streamwise wavenumber demonstrate this same 
structure. A characteristic consistent of this subclass is that the Reynolds 
stress is generated near the wall. The individual velocity components of the 
(1,2,1) mode, the most energetic wall mode, is shown in Figures lllbl and 
llldl for their real and imaginary components, respectively, and in Figures lllcl 
and lllal for their amplitude squared and phase time history. The amplitude 
history shows the bursting nature of these travelling-waves, and the near con- 
stant phase velocity shows why these modes are referred to as propagating or 
travelling-waves. These wall functions constitute 35.22% of the total energy of 
the flow, and are the most energetic of the propagating modes. 

The next mode of interest, still in the propagating mode class, is the lift 
mode, found whenever the streamwise wavenumber is greater than or equal to 
the spanwise wavenumber and the spanwise wavenumber is greater than one. 
Three examples of this subclass are shown in Figures I12al Il2fl and |12g| for 



modes (2,2,1), (3,2,1), and (3,3,1), respectively. In a lift mode, the coherent 
vorticity and Reynolds stress extends from the wall to near the centreline, and 
results in a lifting motion. The real and imaginary velocity components are 
shown in Figures fl2bl and Il2dl respectively, and the amplitude and phase time 
history in Figures I12cl and I12e^ respectively. Again, the constant phase speed 
classifies these structures as propagating modes, and the amplitude bursts 
are also characteristic of a turbulent flow field. The lift structures constitute 
29.68% of the total energy of the flow. 

The next modes, the asymmetric modes, are responsible for breaking the 
symmetry of the flow about the axis of the pipe. These modes are found for a 
spanwise wavenumber of one, for any streamwise wavenumber, and consist of a 
coherent vortex just outside the log layer. They are asymmetric because of their 
nonzero radial and azimuthal velocities at the origin, which is physical only 
for azimuthal wavenumber n = 1, since a positive radial velocity at the origin, 
under a rotation of n around the axis, results in a negative radial velocity. The 
same is true for the azimuthal velocity at the origin. Three examples of these 
modes are shown in Figures Il3al Il3fl and |13g| for modes (1,1,1), (2,1,1), and 
(3,1,1), respectively, again having the Reynolds stress between the coherent 
vorticity. These modes are also propagating and turbulent, as seen in the 
(1,1,1) mode's amplitude (Figure [13cp and phase time history (Figure [13ep . 
with its real and imaginary component found in Figure I13al and I13al This 
mode constitutes 9.09% of the total energy of the flow. 

The last of the propagating modes is the ring mode, named for the ring of 
coherent vorticity that is found for all modes with zero azimuthal wavenum- 
ber. Examples of this structure are shown in Figures I14al Il4fl and 14g for 
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modes (1,0,1), (1,0,2), and (2,0,1), respectively. Since the ring modes have 
no azimuthal dependance, only a radial cross-section of coherent vorticity is 
shown with Reynolds stress superimposed. The real and imaginary velocities 
are shown in Figures Il4bl and Il4d[ and the amplitude and phase in I14al and 
I14el This mode constitutes 6.60% of the total energy of the flow. 

The last two modes are the non propagating modes. The first and more 
energetic of the two is the roll mode, found for any mode with zero stream- 
wise wavenumber. This consists of rolls of coherent vorticity, as seen in the 
three example modes (0,6,1), (0,5,1), and (0,3,1) in Figures fl5al I15fl and|15g 



respectively. The Reynolds stress is strong between the coherent vortices but 
is an order of magnitude less than those of the wall or lift structures. The 
velocity components of the (0,6,1) mode are shown in Figure [T5bl and fl~5dl Of 
note, these structures do not have a constant phase velocity, seen in Figure 
I15et and the decay rate of the energy is slower than that of the propagating 
modes, seen if Figure ll5cl This subclass constitutes 18.34% of the total energy 
of the flow. 

The other non propagating mode is the shear mode, found for zero azimuthal 
and streamwise wavenumber, therefore corresponding to the fluctuation of the 
mean flow rate. Since these structures have no coherent vorticity nor imaginary 
components, only the real velocity components (Figures I16al and ll6cp and their 
amplitude and phase time history (Figures Il6bl and Il6dj) are plotted. Like the 
roll modes, the shear modes do not have a constant phase velocity, and since 
they do not have an imaginary component, the phase oscillates between zero 
and 7r. The shear modes constitute 1.08% of the total energy of the flow. 

The energy spectra of the propagating mode subclasses are shown in Figure 
[T71 This shows that the tail wavenumber end of the inertial range is popu- 
lated primarily by the lift modes and the low end of the spectra is populated 
predominantly by the wall modes. The physical meaning of this is that the 
energy starts at the wall with large scale structures and is lifted away from 
the wall into the outer region in small scale structures, represented by the 
lifting modes. 

The effect of higher radial quantum number q > 1 is more zero crossings 
of the velocities in the radial direction. The effect on the location of coherent 
vorticity and the number of vortex cores scales with the radial quantum num- 
ber q, making the subclasses invariant with q. The wall and lift modes retain 
their characteristics, in that even with more vortices, they remain close to the 
wall for the wall mode and close to the centreline for the lift modes, shown in 
Figures I18al through I19b[ 
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Figure 11a. Most energetic propagating wall mode (1,5,1). Coherent vorticity (left) and a 
cross-section of coherent vorticity with Reynolds stress superimposed (right). 
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Figure lib. Real component of (1,5,1). 




Figure 11c. Time history amplitude of (1,5,1). 
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Figure lid. Imaginary component of (1,5,1). 



Figure lie. Time history phase of (1,5,1). 
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Figure llf. Propagating wall mode (1,3,1). 
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Figure llg. Propagating wall mode (2,4,1). 
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Figure 12a. Most energetic propagating lift mode (2,2,1). Coherent vorticity (left) and a 
cross-section of coherent vorticity with Reynolds stress superimposed (right). 
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Figure 12b. Real component of (2,2,1) 



Figure 12c. Time history amplitude of (2,2,1) 
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Figure 12d. Imaginary component of (2,2,1) 



Figure 12e. Time history phase of (2,2,1) 
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Figure 12f. Propagating lift mode (3,2,1). 
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Figure 12g. Propagating lift mode (3,3,1). 
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Figure 13a. Most energetic propagating asymmetric mode (1,1,1). Coherent vorticity (left) and a 
cross-section of coherent vorticity with Reynolds stress superimposed (right). 
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Figure 13b. Real component of (1,1,1). Figure 13c. Time history amplitude of (1,1,1). 

0.06 | 1 1 I 




2000 4000 6000 8000 10000 12000 14000 16000 18000 



Figure 13e. Time history phase of (1,1,1). 

Figure 13d. Imaginary component of (1,1,1). 
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Figure 13f. Propagating asymmetric mode (2,1,1). 
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Figure 13g. Propagating asymmetric mode (3,1,1). 
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Figure 14a. Most energetic propagating ring mode (1,0,1) cross-section. Coherent vorticity with 

Reynolds stress superimposed. 
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Figure 14b. Real component of (1,0,1). 



Figure 14c. Time history amplitude of (1,0,1). 
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Figure 14d. Imaginary component of (1,0,1). Figure 14e. Time history phase of (1,0,1). 
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Figure 14f. Propagating ring mode (1,0,2). 
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Figure 14g. Propagating ring mode (2,0,1). 
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Figure 15a. Most energetic non propagating roll mode (0,6,1). Coherent vorticity (left) and a 
cross-section of coherent vorticity with Reynolds stress superimposed (right). 
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Figure 15b. Real component of (0,6,1). Figure 15c. Time history amplitude of (0,6,1). 
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Figure 15d. Imaginary component of (0,6,1). Figure 15e. Time history phase of (0,6,1). 




Figure 15f. Non propagating roll mode (0,5,1). 
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Figure 15g. Non propagating roll mode (0,3,1). 
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Figure 16a. Real component of (0,0,1). 



Figure 16b. Time history amplitude of (0,0,1). 
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Figure 16c. Real component of (0,0,2). 



Figure 16d. Time history phase of (0,0,1). 



4 Discussion 



Although the Karhunen-Loeve method yields a preferred or natural basis for 
turbulence, one must be careful in the conclusions drawn from the results, 
as any structure or feature can be reconstructed from any given orthogonal 
basis. The strength of this method is that it is the most optimal basis, and by 
focusing on the most energetic modes a low order dynamical representation of 
the flow is realised. This low order dynamical system reveals three important 
results, and paints a picture of the energy dynamics in turbulent pipe flow. 

The first result is the constant phase speed of the propagating modes. This 
was also found in studies of turbulent channel flow, as the structures repre- 
sented by the basis advect with a constant group velocity, the same average 
velocity of burst events [13-15]. The normal speed locus of the propagating 
waves is shown in Figure [20] for the propagating modes found in the top 50 
most energetic modes. For this, the phase speed w/||k|| is plotted in the direc- 
tion k/||k||. The locus is nearly circular, which is evidence that these structures 
propagate as a wave packet or envelope that travels with speed of 8.41 U T , the 
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Streamwise Wavenumber 



Figure 17. Energy spectra for the propagating mode subclasses, wall (solid), lift (dashed), 
asymmetric (dash-dot), and rings (dots). The wall modes correspond to the low wavenumber energy 
and the lift modes represent the high wavenumber spectra commonly known as the inertial range. 

point at which the circle intersects the x-axis in Figure 1201 which corresponds 
to the mean velocity in the buffer region at y + = 9.6. 

The second interesting result is that none of the modes with azimuthal wave 
number n = 1 exhibit any travelling waves near the wall, and are without 
exception a streamwise or inclined streamwise vortex in the outer region. The 
reason is that the n = 1 mode does not allow for a near-wall travelling wave 
as found by Kerswell [34], and as such, the basis expansions for these modes 
have only near-centreline structures. 

The third result of note is that the different KL structures, as observed 
from our results, can be seen qualitatively as an expansion that represents the 
horseshoe (hairpin) vortex representation of turbulent structures. This horse- 
shoe structure is supported by a large number of researchers in the turbulence 
community as the self-sustaining mechanism for turbulence [1]. In the repre- 
sentative horseshoe structure (see, for example, the figure by Theodorsen [35]), 
the structures found in the KL can be seen. The wall modes represent the leg 
structure and its perturbation near the wall. The lift modes represent the 
structure lifting off the wall. The asymmetric modes represent the secondary 
and tertiary horseshoes that are formed. The turn modes represent the span- 
wise head of the horseshoe. As mentioned before, since the KL decomposition 
forms a basis, any flow can be recreated from these eigenmodes, so this result 
is reported as an interesting qualitative result, and as an understanding that 
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Figure 18a. (6,2,3) lift mode, contours of 
coherent vorticity. 
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Figure 18b. (6,2,5) lift mode, contours of 
coherent vorticity. 
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Figure 19a. (2,6,3) wall mode, contours of 
coherent vorticity. 



Figure 19b. (2,6,5) wall mode, contours of 
coherent vorticity. 



the propagating modes are the structures of interest in turbulence. 

The final dynamical picture of the KL modes in summarising the work done 
by Kerswell [34], Webber et al. [20,21], Sirovich et al. [14], and this current 
work is as follows. Webber et al. found that the energy enters the flow from 
the pressure gradient to the shear modes. The shear modes, interacting with 
the roll modes present in fully developed turbulence transfers the energy from 
the shear modes to the roll modes. Shown in this study, and theoretically by 
Kerswell, the roll modes decay much more slowly than the propagating modes. 
Because of this slow decay, these streamwise rolls provide an energy storage 
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Figure 20. Normal speed locus of top 25 most energetic propagating eigenfunctions, normalised 
with shear velocity u T . The circle is a least-squares fit to the data and represents that the wave 
packets are acting together as a group with speed of 8.41 U T , the point where the circle intersects 

the y-axis. 



role in the turbulence, similar to a capacitance. This allows enough time for 
the roll to interact with the propagating wave packet shown by Sirovich et al. 
Because this interaction requires a finite time to occur, the energy storage and 
slow decaying nature of the streamwise rolls play an integral role in the self- 
sustaining nature of turbulence. It is also surmised that the travelling modes 
found by Kerswell and other researchers are represented in the KL formulation 
as the most unstable wave packet of KL modes. Following again the work of 
Webber et al. and applying the classifications found in this study, the majority 
of the energy from the roll modes are transferred to the wall modes. The wall 
modes then interact with themselves through an asymmetric mode catalyst, 
and with the lift modes through a ring mode catalyst. Physically, this gives 
a picture of wall turbulence energy being transferred near the wall and then 
lifted up to the outer layer. This interaction between the wall and the lift 
modes is what populates the inertial range shown in figure [T71 The energy is 
ultimately dissipated by viscosity, faster for the higher wavenumber modes. 
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5 Conclusions 

We have presented the use of the Karhunen-Loeve expansion method with the 
results of a globally high-order direct numerical simulation of turbulent pipe 
flow. The results reveal the structure of the turbulent pipe flow as propagat- 
ing (80.58% total energy) and non propagating modes (19.42% total energy). 
The propagating modes are characterised by a constant phase speed and have 
four distinct classes: wall, lift, asymmetric, and ring modes. These propagat- 
ing modes form a travelling wave envelope, forming a circular, normal-speed 
locus, with advection speed of 8A1U T corresponding to the mean velocity in 
the buffer region near y + = 9.6. The non propagating modes have two sub- 
classes: streamwise roll and shear modes. These represent the energy storage 
and mean fluctuations of the turbulent flow, respectively. The energy is trans- 
ferred from the streamwise rolls to the wall modes first, and then later to 
the lift modes, physically representing the energy transfer from the wall to 
the outer region. This eigenfunction expansion, using both their structure and 
their time-dependent coefficients, provides a framework for understanding the 
dynamics of turbulent pipe flow, and will provide a basis for further analysis 
and comparison leading to understanding the mechanism of drag reduction by 
spanwise wall oscillation [36] and the mechanism of relaminarization [37] . 
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